Simulating Stochastic Dynamics Using Large Time Steps 
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We present a novel approach to investigate the long-time stochastic dynamics of multi-dimensional 
classical systems, in contact with a heat-bath. When the potential energy landscape is rugged, the 
kinetics displays a decoupling of short and long time scales and both Molecular Dynamics (MD) or 
Monte Carlo (MC) simulations are generally inefficient. Using a field theoretic approach, we perform 
analytically the average over the short-time stochastic fluctuations. This way, we obtain an effective 
theory, which generates the same long-time dynamics of the original theory, but has a lower time 
resolution power. Such an approach is used to develop an improved version of the MC algorithm, 
which is particularly suitable to investigate the dynamics of rare conformational transitions. In the 
specific case of molecular systems at room temperature, we show that elementary integration time 
steps used to simulate the effective theory can be chosen a factor ~ 100 larger than those used in 
the original theory. Our results are illustrated and tested on a simple system, characterized by a 
rugged energy landscape. 

I. INTRODUCTION 

The investigation of a vast class of physical phenomena requires the understanding of the long-time dynamics of 
classical systems, in contact with a heat-bath. Examples include critical dynamics, molecular aggregation, protein 
folding, to name a few. 

The most natural strategy to describe these processes is to integrate numerically the equations of motion, i.e. to 
perform MD simulations. Unfortunately, when the number of degrees of freedom is very large, or in the presence 
of large free energy barriers, MD approaches become extremely costljD} or even impracticable. The problem arises 
because the time scale associated with the system's local conformational changes can be many orders of magnitude 
smaller that the time scales of the dynamics one is interesting in studying. As a result, most of the computational 
time is invested in simulating uninteresting thermal oscillations. 

This situation is exemplified in FigJTJ where we show the stochastic motion of a point particle, interacting with a 
2-dimensional external potential. The solid line was obtained by means of a MD simulation and illustrates how, at 
short time scales, the dynamics of this system is dominated by fast modes associated to thermal diffusion. However, 
when the evolution of the system is described using much lower time resolution power, the effect of such short-time 
thermal fluctuations tends to average out and to become unimportant. This is evident from the comparison between 
the solid line and the dashed line, which was obtained by averaging over blocks of consecutive frames in the original 
MD trajectory. At long times, the dynamics of system is mostly sensitive to the structure of the external energy 
landscape, which was chosen to be spherically symmetric. 

Clearly, an important question to ask is whether it is possible to develop theoretical/computational frameworks 
which yield directly the correct long-dynamics, but avoid investing computational time in simulating the short-time 
thermal oscillations. Significant progress in this direction has been recently made developing approaches based on 
Markov State Models^EEl A potential difficulty of such approaches resides in correct identification of the metastable 
states. In addition, for each different system, one needs to perform a large set of independent MD simulations in order 
to accurate calculation of the rate coefficients. 

In this work, we present an alternative approach to simulate the dynamics over long times. We develop a rigorous 
effective theory which (i) yields by construction the correct long-dynamics and (ii) does not require to identify meta- 
stable states, nor to evaluate the transition matrix by MD. To our goal, we use a field theory approach, based on 
Renormalization Group (RG) ideas and on the notion of effective field theorjEl Such a powerful tools have been 
already successfully applied to describe the low-energy dy namics of a vast variety of of quantum and statistical 
systems characterized by a separation of scales — see e.gP^ — . To the best of our knowledge, this method has never 
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FIG. 1: Langevin dynamics of a point particle in a 2-dimensional external potential. The solid line denotes the result of an MD 
simulation. The dashed line is the result of averaging over blocks of consecutive frames of the MD trajectory. Such an average 
smoothes out the trajectory. 

been applied to develop an effective theory to efficiently simulate the long-time stochastic dynamics of a system in 
contact with a heat bath. 

The main idea of our approach is to exploit the decoupling of time scales in the system in order to define a pertur- 
bative series, in which the expansion parameter is the ratio of short- over large- time scales. In such a perturbative 
framework, the average over the short-time fluctuations can be computed analytically, to any desired level of accuracy. 
The average over the fast thermal oscillations gives rise to new terms in the stochastic path integral, which represent 
corrections both to the interaction and to the diffusion coefficient. Such new terms implicitly take into account of the 
dynamics of the fast degrees of freedom, which have been integrated out from the system. 

Once a finite number of such effective terms corresponding to a given accuracy have been calculated analytically, 
it is possible to simulate the dynamics of the system using much larger time steps. By construction, in the regime of 
decoupling of fast and slow modes, one is guaranteed that the effective long time theory generates the same probability 
distributions of the underlying, more fundamental stochastic theory. It is important to emphasize the fact that the 
present approach is not equivalent to simply including higher-order corrections in the Trotter expansion 8 ! Indeed, the 
assumption of decoupling of time scales leads to further simplifications with respect to such an approach. 

The paper is organized as follows. In section [IT] we review the path integral formulation of the Langevin dynamics 
and we outline the formal connection between stochastic dynamics and evolution of a quantum particle in imaginary 
time. Such a connection is used in section [TTT] to identify and isolate the dynamics of the fast degrees of freedom. In 
sections |IV| and [V] we present our perturbative scheme which allows to integrate out the fast modes and derive the 
effective interactions and diffusion coefficients. In section |VI| we discuss how the effective theory for the dynamics 
at long time scales can be simulated using the diffusion MC algorithm, which is briefly reviewed in appendix [5] 
Section |Vl| is devoted to simple examples, which illustrate how this method works in practice. In section |YlII| we 
discuss the applicability of the present approach to simulate the Langevin dynamics of molecular systems. Results 
and conclusions are summarized in section llXl 



II. LANGEVIN DYNAMICS 

We consider a system defined by a stochastic d-dimensional variable x obeying the Langevin Eq.n: 

mi = -VU(x) ~ 7± + £(t), 



(1) 
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where U(x) is a potential energy function, m is the mass, 7 is the friction coefficient and is a S— correlated 
Gaussian noise. In many molecular systems of interest, the acceleration term mx is damped at time scales of the 
order 10 -13 s, which much smaller than the time scale associated to local conformational changes. If such a term is 
dropped one obtains the so-called over-damped or velocity Langevin Eq.: 

x = --VU(x)+r)(t), (2) 
7 

where rj(t) is a rescaled delta-correlated Gaussian noise, satisfying the fluctuation-dissipation relationship: 

( V ( t ')r ] (t)) = 2dj^5(t-t'). (3) 

The over-damped Langevin Eq. defines a Markovian process. The probability distribution P(x,t) generated by 
such a stochastic differential equation obeys the Fokker-Planck Eq. 

%-P(x, t) = i V [VP(x, t) + pVU(x)P(x, t)] . (4) 
at p 7 

By performing the substitution 

P{x,t) = e~i u{x) fl>(x,t) (5) 
the Fokker-Planck Eq.Q can be recast in the form of a Schrodinger Equation in imaginary time: 

- ~V>0M) =H eff ip(x,t), (6) 
where the effective "Quantum Hamiltonian" operator reads 

Heff = -~V 2 +(3V eff (x), (7) 
and V e ff(x) is called the effective potential and reads 

Veff(x) = i ({VU(x)f - -^ 2 U{x)^ . (8) 

Hence, the problem of studying the diffusion of a classical particle can be mapped into the problem of determining the 
quantum-mechanical propagation in imaginary time of a virtual system, defined by the effective quantum Hamiltonian 
([7]), interacting with the effective potential V e ff(x). 
Let G(xf,tf\xi) be the Green's function of the Fokker-Planck operator, subject to initial condition x(0) = a;,-, i.e. 

~G(xf,t\xi) - j- V [VG(x f ,t\ Xl ) + pVUG(xt,t\xi)] = S(t)S(x - Xi ) (9) 

The interpretation of such a Green's function is the probability for the system to be in x at i, conditioned to start 
from Xi at the initial time. Formally, such a conditional probability can be related to the "quantum" propagator of 
the effective Hamiltonian (|7j: 



G{x,t\xi) = exp 



l(U{x)-U{ Xi )) 



K(x,t\ Xi ), (10) 



K{x,t\xi) = (x\e- tH °ff\xi). (11) 
Hence, it is immediate to derive a path integral representation of the Green's function G(x, t\xi): 

G(x,t\ Xi ) = e -f>/W(*)-u(*,) j Xit) ~ X Vxe-P s 'ffW, (12) 

J x{ti)—Xi 



where S e f / [x] is the effective " action" , 



S eff [x]= f\r Qi 2 +V eff (xj). (13) 
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The pre-factor e ' 3 / 2 C 7 ( a: ) u i x t) in Eq. |llj can be transformed away, noticing that = xll'(x). One than obtains 

a path integral in which the statistical weight contains the Onsager-Machlup functional 



G(x,t\xi) 



x(t)—x 



Vx e 



-fS f'dr U x 2 +i x U'(x)+V ef f(x)j. 



x(0)=Xj 



(14) 



Eq. ( 12 ) provides an expression for the conditional probability in terms of the microscopic stoc hastic dy namics 
governing the system. It represents the starting point of the Dominant Reaction Pathway approac h 1 9 * 10 * 11 * 12 -* which 
deals with the problem of finding the most probable transition pathways between the given configurations Xi and Xf, 
which are visited at the initial and final time x(t) = Xf, x(0) = Xi, respectively. 

On the other hand, in this work we are interested in the corresponding initial value problem, i.e. we are want to 
develop an effective theory which yields directly the long-time evolution of the probability density P(x,t), solution 
of Eq. Q, starting from a given initial probability density P(x,t = 0) = Pq(x). The probability density P(x,t), the 
Green's function G(xf,t\xi,ti) and the initial distribution po(x) are related by the Eq. 

P(x,t) = [dy G(x,t\y) p (y). 



(15) 

Hence, for positive time intervals, the conditional probability G(x, t\xi) can be considered as the propagator associated 
to the stochastic Fokker-Planck Eq. Q. 



III. SEPARATION OF FAST AND SLOW MODES 



Without loss of generality, let us focus on the stochastic path integral (12 1, with periodic boundary conditions: 



Z(t) = / dx G(x,t\x,0) 



T>x exp 



(16) 



We observe that the inverse temperature ^ plays the role of H, in the analogy with the quantum mechanical formalism. 
Hence, the loop expansion of the path integral (16 1 generates an expansion in powers of jj. 
Let us introduce the Fourier conjugate: 



x(w n ) 



1 /■' 

- / drexp [— iui n t) x(t) 
1 io 



j(t) = x(t + t) = x{uj n ) exp [iu! n t] 



(17) 
(18) 



where uj n are the Matsubara frequencies: 



2tt 

T 



n = 0,±1,±2, 



(19) 



In numerical simulations, the integration of the over-damped Langcvin Eq. is performed by choosing a finite 
elementary time step At. In frequency space, this implies the existence of an ultra-violet cut-off £1, which is inversely 
proportional to At: 



9. 



2tt 
At' 



(20) 



Such a relationship becomes a strict equality in the case of periodic boundary conditions, as in Eq. (16 1. In general, 



when the boundary conditions are not periodic, it represents just an order-of-magnitude estimate of the largest Fourier 
frequencies, which are associated to a given choice of the integration time step At. 

Let us now introduce a parameter < b < 1 and split the frequency interval (0,f2) as (0,6 O) U (b 0,0). Then 
the Fourier decomposition of a path contributing to (IT6|) can be split as: 



x(t) — x > {t) + x < (t), 



(21) 
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where x < (t) and x > {t) will be referred to as the slow- and fast- modes respectively: 

x K (t) = J2 i K) e '""' 

\uj n \<bQ 

x> (t) = Yl %) ei " n< 

bn<\u)„\<Q 



(22) 
(23) 



The main purpose of this work is to develop a perturbation scries to systematically integrate out from the path 
integral the modes with frequencies uj n > ML To this end, we begin by re-writing the "kinetic" term which appears 
in the effective action ( 13 1 of the path integral ( 12 1 as a sum of the kinetic energy of slow and fast modes: 



dr x 



tjy 
4 

7 
4 



^ x(uj n ) X (-L) n ) 

|w„|<0 



dr x 2 Jt) + 



It 



^2 w « 



(24) 



where Sb denotes the shell of hard modes Sb — (M7, ft). 

Let us now consider the potential term and expand around the slow modes £<(£) 



V eff (x(r)) = V eU {x < {r)) + 



dV eff ( X< (r)) 
dx l 



1 d 2 V eff ( X< (T)) 

2 dx l x : > 



4(r) xi(r)+0(xi) 



(25) 



The complete path integral (16 1 can be split in the following way: 

Z(t) = j>Vx < j>Vx > e -P s eff [ x< (t)+ x> (t)] 

= j>V X< e-/ 3 ^// [-<(*)] e -^>[x<(r)]_ (26) 

in this expression, the action functional 5* e // is evaluated on the slow-modes only and depends on the original effective 
potential V e ff (which we also shall refer to as the "tree-level" effective potential). S > [x < (t)] is a correction term 
action which accounts for the dynamics of the fast modes which are integrated out: 



-/3S>[x<(r)] 



V X> e 



i< 5(w„) x(-u> n )-0S irl t 



(27) 



where the Si n t is an effective interaction term. In such an Eq., the integration over the hard modes is performed in 
Fourier space, 



Vx > = ] [ dx(uj n ). 
|«„| es b 



(28) 



Eq. (26 1 is formally exact. In the next section, we evaluate the effective action S > [x < (t)] perturbatively. The 
effective interaction which includes the correction coming from S' > [x < ] will be referred to as the renormalized effective 
interaction. 



IV. RENORMALIZED EFFECTIVE INTERACTION 



In the previous section, we have seen that the the integration over the fast modes generates an additional term in 
the effective action for the slow modes: 

Z(i)s ^ <e -,,,H.<<,,,-^< W1 (29) 

where 

e - s>[x<(t)] = I Vx> e -/« E,„ nleS6 ^ *(— ») e -ps int (30) 



(1PR): 




V1 




V3 



(1PI-D): 




Sb 



(1PI-nD): (v2) 





V2 



V4 



FIG. 2: Examples of connected graphs appearing in the exponent of Eq. (38 1. The diagrams on the upper part (dumbbell 



diagrams) are one-particle reducible, while those in the middle and in the bottom are one-particle-irreducible. In particular, 
those in the middle (daisy diagram) are local in time. 



In this section we formally perform such an integration. We begin by re-writing e & s >l x <( T )} as 

e -0 S>[*<(r)] = ( e -^"-) , 

where the notation (-)o denotes the expectation value evaluated in the free theory 

,2 



t ^2 %( w n) x(-u n ). 



To evaluate the matrix clement (e /3Si " t ), we represent the e P Sint "operator" by its power series: 

k > k 

Next, we expand the interaction action Si n t[x > + x<] around the slow modes 1 

'•' dv eff ( x< ( T )) 



i r , d 2 v eff ( x< (T)) 



ps mt [x > + x< ] = -pj dr " 11 ,■ ( 7" i — / ,/- " in/-' (7-) + .. 



dx l 



where i h (r) £>(t) . . . x> (r) are vertices with couplings 



Vi u ...,i h (r) 



a fc K//[x<(r)] ^ 
dx 11 . . . dx lk 



(31) 



(32) 



(33) 



(34) 



(35) 



Throughout all this work, we shall adopt Einstein notation, i.e. the summation over repeated indexes is implicitly assumed. 
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FIG. 3: Diagrammatic representation of the local time-derivative expansion of a non-local diagram — Eq. (50 1 — . Solid lines 



are fast-mode propagators, while dashed lines represent a single time derivative acting on the corresponding vertex function. 



Notice that each term in the perturbative expansion ( 34 1 generates a new vertex, with an increasing power of the 



x > (t) field. The couplings to the fast modes depend implicitly on the time r, through the slow modes x < (r). 



By Wick theorem, each term in the series (33 1 can be related to a Feynman graph with vertexes given by (35 1 and 
propagators given by — see appendix [A] — : 



(xKn) a^(r 2 ))o 



\u m \,\u n \£Sb 



Lr > (U> n ,U} m ) e — 



l(T2-Tl) 



\Un\£bb 



(36) 



The expansion (33 1 can be re-organized as the exponent of the sum performed over only connected diagrams: 



/3 S>[a:<(r)] ^(sum over all connected diagrams) 



Hence, the path integral (26) for the slow modes can be given the following exact diagrammatic representation 

-/? S , e //[a;<(i)]] + (sum over all connected diagrams) 



z{t) 



T>x < e 



(37) 



(38) 



Below we give a classification of all the connected diagrams that may give a contribution to the expansion above. 
Firstly note that all diagrams that involve an odd numbers of fast field vanish thanks to the Wick theorem. We are 
thus left with the following sets of, a priori nonvanishing, diagrams: 



• 1PR (one-particle-reducible) diagrams, namely diagrams that can be topologically separated into two distinct 
subdiagrams by cutting one internal fast mode line (propagator): they have the topology of a dumbbell. The 
simplest examples of dumbbell diagrams are depicted in the upper part of Fig. [2| 

The main assumption of this work is the existence of a gap between slow modes and fast modes. Under such 
assumption all the 1PR diagrams give vanishing contributions. From the physical point of view, this can be 
understood as a consequence of energy conservation: in order for the total energy flowing through a vertex with 
a single hard mode to be conserved, at least one of the external modes has to be hard. On the other hand, 
our working assumption implies that all the modes in the external legs of diagrams are soft. This result can be 
rigorously proven for all 1PR. As an example, we explicitly compute the upper left diagram of Fig. [2j We have 



1 / dry \ dr 2 V t (x<(n)) Vj (x<(r 2 )) £ — 



2! V 2! 



(39) 



We note that the effective potentials depend smoothly on time, through the periodic functions x<(r). Hence, 
the terms Vi[x<(n)] and Vj[x<(r2)] in Eq.(39l can be expressed in terms of their Fourier-transform, 



^Or<(n)) 



n 



(40) 
(41) 



This allows to perform the time integrals, which simply yield t 2 S UJn+ ^ n ,o S Un - Vm ja. Due to such delta-functions, 
only hard z^-modes survives, which are projected in a term 



E 



1 



fit 
4 7 »l »l 



Vi(v n )Vi(-V n ) W 0. 



(42) 
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These modes thus yield negligible contributions under the physical assumption of large separation of frequency 
scales. On the other hand, if one does not assume a separation of time scale, this diagram gives finite contribution 
and has to be accounted for. Note that this term has the same structure as the first correction which appears 
when one performs higher-order Trotter expansiorP. 

It is not difficult to check that such result holds for all 1PR diagrams, so that we can reduce our effective action 
to the sum of 1PI (one-particle-irreducible) diagrams, i.e. diagrams that cannot be disconnected by cutting a 
single internal line. They can be classified in two main groups: 

• 1PI "daisy" diagrams, namely diagrams with a single vertex. Such diagrams only involve equal-time hard 
propagators and only give rise to contributions to the renormalized effective action which are local in time: they 
have the topology of a daisy, hence the name. Examples of daisy diagrams are depicted in the middle part of 
Fig.[2j It is not difficult to compute a generic daisy diagram with K petals (propagators). It is due to the vertex 
with 2K hard fields and reads 



(43) 



where A = <£y didj is the Laplacian operator, and the numerical factor in front is a combinatorial factor. The 
sum, i.e. the equal-time propagator, can be easily performed by taking the continuum limit — ► ^ J duj that 
simply yields 2 

1 y J^ 1 / Qrf " = 1 x - h2p - x (44) 



so that we finally obtain 



0\ (Dl-b 



K r t 



o 



dr A^y e// (x<(T)), (45) 



where we have reinstated the diffusion coefficient D = 1/(3"/. Hence, one can even formally resum all the daisy 
diagrams into the compact expression 



^ daisy diagrams = —(3 J dr exp ( ^^-fS 



(x<(r)) . (46) 



• 1PI non-daisy diagrams: all other non-local diagrams. The simplest examples of such diagrams are depicted in 
the lower part of Fig. [2] These diagrams generate contributions to the renormalized effective action that are 
non-local in time and give rise to infinite series of local diagrams. For example, the evaluation of the lower left 
diagram of Fig. [2] yields a contribution of the form: 

2 x 2- (-1) I d ^ I dT * v - ( *< (Ti)) Vki ( *< (T2)) £ {^t ) — ^ 6 



'Hk Vjfo 



n 



(47) 

where the 2 in front is a combinatorial factor. After Fourier transforming the potentials (see discussion below 



eq. (39)), the integrals over times yield t ^ m +w„+i/ 7l ,o <5u m +w„-i> m .o- Hence, 



7?,?. i *'^' J 3k i tf' (48) 



|^„|es b 



2 



Here for later use we consider a generic even power 2p. It is easy to check that the error one makes in considering the continuum limit 
is of order nnMlv with N = Qt/2ir. 
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Now we again make use of the assumption that slow modes and fast modes of physical processes under study are 
separated by a large gap. Under such assumption we can safely expand the second fraction in the latter expression 
in power series of slow modes v n and rewrite (48) as higher-time-derivative expansion. Let us reintroduce the 
integral over time as 1 = \ J* dr e l (»ri+u m )T so that powers of v n can be traded with time-derivative of the 

potential (note that odd powers vanish upon symmetric sum; in fact they would give rise to total time-derivative 
terms that are zero upon integration thanks to periodicity in time.) We are thus left with 



2 / dT ^0<( T )) 

7 Jo 



L \u n \es b 



|w„|es 6 



\u n \es b 



(49) 



Sums over hard frequencies can be performed in the continuum limit with the help of formula (44 1 and time- 
derivatives can be partially integrated in order to rewrite the latter in a more symmetric form 



1 



Ti" V Jo 



dT 



1 1 - b 3 

3 (&fi) 3 



(Vij(x<) 



2 3 1 - b 5 



2 5 1 - b 7 



(50) 



The infinite higher-derivative expansion is the legacy of non-locality in time: such an expansion can diagram- 
matically represented as an infinite sum of local (daisy-like) diagrams, as depicted in Fig. [3j 

It is intuitive to expect that, in the presence of decoupling low and high frequency modes, the higher-derivative terms 
should be suppressed. In the next section, we shall generalize this statement and present a quantitative method to 
systematically organize all contributions to the effective action in terms of a perturbative series. 

V. SLOW-MODE PERTURBATION THEORY 



The diagrammatic representation of the path integral given by Eq. (38 1 is formally exact, but rather useless. In 



fact, it is obviously impossible to evaluate and re-sum exactly all the infinitely many Fcynmann graphs appearing in 
the exponent. On the other hand, in this section we show that it is possible to compute the renormalized effective 
action S e f f[x < (t)]] to an arbitrary level of precision, by calculating only a, finite number of Feynmann graphs. This 
way, the low-frequency effective theory retains predictive power. 

The idea is to exploit the decoupling of the short-time dynamics from the long-time dynamics to organize the sum 
over all possible graphs as a perturbative expansion. We shall refer to such a systematic evaluation of the renormalized 
low-frequency effective action as to the slow-mode perturbation theory. 

The first step in the construction of our perturbation series is to identify all the dimcnsionless combinations of 
the physical quantities which appear in the Feynmann graphs contributing to ( 38 1 , evaluated in stationary phase 
approximation. Let us first define the quantities 



V 



dr V 



V, 



k 



1 

f 



dr A m V 



1 

U 2 



dT 



(d™vy 



(51) 



where k is the typical wave vector on the spatial Fourier transform of V e ff{x) and lu < is the typical frequency in 
temporal Fourier transform of V e ff(x < (r)). 

Using these combinations, we can thus construct the following dimensionless combinations: 



0V 

bn 



a-2 



k 2 D 



a 3 



b9. 



(52) 



We are interested in describing the dynamics of physical systems for which each of these parameters can be considered 
small. In order to illustrate the physical interpretation of the condition a\ <C 1, we observe that the probability for 
the system to remain in the same configuration x, during an elementary time interval dt is 



P(x, t + dt\x, t) oc 



1 



(dty/ 2 



-0V eff (x)dt 



(53) 



Hence, the combination (3V e ff represents 3 the typical time scale associated to local conformational changes, and the 



3 Notice that, in the small temperature limit, V e ff(x) becomes positive definite. Thus, P(x, t + dt\x, t) decays exponentially with PV e ff 
in the time interval dt. 
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(L=2) 



ceo cego 



(L=2) 
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FIG. 4: Examples of the diagrams with the lowest degree of slowness, up to L — 2. 

condition a\ <C 1 expresses the condition that the time spent on average by the system in each configuration is large 
compared to the elementary short-time scale, dt ~ -r^. 

The condition <C 1 implies that the effective potential varies over length scales which are large, compared with 
the mean distance covered by Brownian motion in an elementary time interval dt. Finally, the condition 03 < 1 
implies that the typical slow mode frequencies are small compared to the ultra-violet cut-off, which is of the order of 
the inverse of the elementary time interval dt. 

It is easy to see that any local diagram in the expansion of the renormalized effective action comes about with 
integer powers of these coefficients, when compared to the tree level effective action. In particular, any diagram 
composed by r vertices of Mi, . . . ,M r hard fields will involve M = XOi spatial derivatives and 4f propagators 
each of which yields a power of . Finally, each additional vertex yields a power of f3V and each time derivative 

_i — 

yields a power of uj. So, the above diagram, at the lowest level in time derivatives will be of order a\ a 2 2 with 
respect to the tree level expression. Higher time derivative terms will add powers 03. It is thus natural to define a 
degree of slowness L for a local diagram, given by 

i(Feynman diagram) = N v — 1 + N T H — (54) 
where A/^ is the number of vertices, N T the number of time derivatives and N x the number of spatial derivatives. The 



definition in Eq. (54 1 is normalized in such a way that L(tree) = 0. 

It is easy to check that the degree of slowness L corresponds to the power of ^ of the local diagrams. Note also 
that for daisy diagrams and for all other diagrams where N T = 0, the degree L is nothing but the number of loops. 
One can thus easily write down and compute the finite set of local diagrams that renormalize the effective action up 
to a fixed (yet arbitrary) level of precision L max . Let us consider a few simple examples. 

L < 1 corresponds to a single daisy diagram with L = 1, jV"„ = 1 and N x — 2, represented in the left panel of 



Fig|4] The expression of this diagram is given by Eq. (45 1 and gives a correction to the effective action of the 
form 



,s [,■ :L, = D J\J J q dr AK ; //(x<(r)). (55) 



L < 2 corresponds to two further diagrams, one daisy diagrams with either N x — 4 and the two-vertex local 
diagram with N x — 2 and no time-derivatives. This latter however is 1PR and gives no contribution. We are 
thus left with the corrections 

2] = ^^0 J\ T Ay e// (x<(r)) + i J* dr *?V eJ j{x k {t)) (56) 
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• L < 3 corresponds to two further diagrams, one daisy diagrams with N x = 6 and the two- vertex local diagram 
with N x — 4 and no time-derivatives. This latter can be simply read off from eq. (50 1. Hence 



S>[x<;L< 3] 



D 1 



7T bfl 

1 /Dl-b 
3! i W 



dr AV eff (x < (r)) + 



1 fDl 



2! 



dr A 3 T/ e// (x<(r)) 



7T 6ft 

D 2 1 - 6 3 
3tt (6fi) 3 7 



/ dr A 2 K//^<(r)) 
Jo 

dr (a 4 a j 14 // (x < (r))) 2 (57) 



that involves in the last term the trace of the square of Hessian of the tree level potential (didj V(x < )) 2 — Tr Tiy. 
In order to see the first time derivatives appearing into the renormalized effective action we need to consider 
L < 5 where, along with several other corrections, we have the correction coming from the second term in (50 1 
that yields 



3/3L> 2 1 - 6 5 
5V (6ft) 5 J 

that can be also recast as a correction of the kinetic action 



dr Tr U\ 



/?7 / dr 



1 3f3 2 D 3 1 - b\ ta 



(58) 



(59) 



Some comments on the results obtained in this section are in order. First of all, we emphasize that the effective 
interactions have been derived under the assumption that the modes which are relevant for the long-time dynamics 
vary over time scales much longer than that of the fast modes, which enter in the loop diagrams. This is the crucial 
assumption of all renormalization group approaches. Our results confirm the intuitive picture that if one adopts a 
low " time-resolution power" , then the effective interactions generated by the ultra-violet modes can be regarded as 
instantaneous. This is in fact general property of renormalization group theory, which is preserved to any order in the 
perturbative expansion. Finally, we note that the correction terms generated by the integration over the fast modes 
is suppressed, in the small temperature limit. 



VI. RENORMALIZATION GROUP IMPROVED MONTE CARLO 



The usefulness of the renormalization procedure resides in the fact that it gives rise to an effective theory, in which 
the largest frequency scale is lowered form ft to 6ft. Equivalently, the shortest time scale is increased form At to | At. 
By construction, in the regime of decoupling of fast and slow modes, the probability density generated by the new 
slow-mode effective theory must be the same as that of the original (i.e. tree-level) theory. In this section, we show 
how it is possible to use the slow-mode effective theory to develop improved MC algorithms for the time evolution of 
the probability density P(x,t), in which the elementary time step used to propagate the configurations is increased 
by a factor 1/6. 

The starting point of the MC approach 13 is to write the probability of observing the system in configuration x at 
time t in terms of the Green's function of the Fokker-Planck Eq. G(x,t\x i ,t i ): 

P(x,t) — J dxi G(x,t\x l ,t i ) po(x l ) (60) 

where pa(y) is the density of states at the initial time. 

One then uses Trotter's formula to write the transition probability as a sequence of intermediate elementary prop- 
agation steps: 

. N-l 

P(x,t)= Y[ dy k G(y k+ i,t k+ i\y k ,t k ) po(xi) (y = x i: y N = x) (61) 

J k=0 

If a sufficiently large number of intermediate steps N is adopted, then the time steps At = t k+ \ — t k can be considered 
infinitesimal and the (unnormalized) transition probability G{y k j r \, t k+ \\y k , t k ) can be calculated analytically 

G(y + dy : t + At\y,t) = const, x e _/3 (^(^) 2 Ai +2^' vc/ fe)) A * e -P v e //(y)At (62) 
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Ito vs Stratonovic in Diffusion MC 

Average posilion at thermal equilibrium 




FIG. 5: Average position at thermal equilibrium, obtained from diffusion MC simulations with (circles) and without (square) 
the branching factor W(t) of Eq. (66 1, for different values of the discretization time step At. Errors are smaller than the 
symbols. 



"Completing the square" in the first exponent, one finds 

G(y + dy, t + At\y, t) = const, x ( dy+ ^ ^( y )f + ^ufAt -p v eff (v) At (63) 

Now and recalling the definition of the effective potential ^ in the second exponent, this Green's function can be 
written as: 

G(y + dy,t + At\y, t) = const, x ( dy+ ^ vu{y) ¥ e &' uw At (64) 

In the MC algorithm, one starts from a set of initial system's configurations, sampled according to he distribution 
Po(xi). Such an ensemble is evolved in time, according to the following procedure. Each configuration is propagated 
for an elementary time interval At, by sampling from the Gaussian 



in Eq. (64 1. Such a configuration is then re-weighted according to the factor 

W(y) = v2 ^) At . (66) 

The iteration of such a procedure for many consecutive elementary propagations gives rise to a set of diffusive 
trajectories, called walkers. In the so-called diffusion MC algorithm, the term W is used to replicate or annihilate 
the walkers. The ensemble of configurations obtained according to this procedure is distributed according to the 



probability density ( 60 ) 



For the MC algorithm to be efficient, the fluctuations in the statistical weight of the walkers — or, cquivalently, in the 
number of walkers — should remain small, throughout the entire time-evolution. This condition is verified if the factor 
W(y) is always of order one. Note however that this term tends to enhance (suppress) the weight of configurations in 
the vicinity of the local minima (maxima) of U (y) , where the Laplacian is positive (negative) . Hence, if the energy 
landscape varies very rapidly in space, then the fluctuations in the statistical weights — or in the number of walkers — 
will in general be large, unless the elementary time step At is chosen very small. This feature represents a limiting 
factor of MC simulations, which makes the sampling of the probability density at large times very computationally 
expensive. 

Clearly, the elementary propagation time step At is the shortest time scale in the simulation. On the other hand, 
in the slow-mode effective theory one integrates out the dynamics in the time scale range (At, 1/bAt). Hence, we 
expect that by taking into account of the corrections associated to the renormalized effective interaction it is possible 
to perform MC simulations in which the integration time step At is chosen a factor 1/6 larger. 
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In practice, the effective slow-mode theory introduces a correction in the re-weighting — or branching — . To order 
L = 1 one has 

Wi=i(y)= e 2v2(J ^ At x e -f ^ vV c// (y)At (67) 

Notice that this expression contains a factor of the inverse frequency cut-off 1 /SI in the exponent. Such a term is 
proportional to the elementary time step At. The corresponding proportionality factor reads 2n only for periodic 
path integral. For a generic initial value MC one can write in general 

Sl = K~ t , (68) 
where the constant K is to be determined from simulations. Hence, we obtain 

W L=1 (y) = e 2 ^ u ^ At x e -« ^ V v*v cff ( y )At> _ (69) 

The unknown constant k can be determined by matching the results obtained by running a short simulation in 



the tree-level theory — i.e. using an integration step At and the tree-level weighting term (66 1 — with those obtained 
in the effective theory — i.e. using an integration step 1/6 At and the renormalized weighting term (69 1 — . In the 
regime of decoupling of fast and slow modes, once the matching has been done, the two algorithms must generate the 
same evolution for the probability density at any later times. 

In the next session, we shall provide an example which illustrates how this procedure works in practice and show 
that the fundamental and the effective theory do indeed generate the same long-time stochastic dynamics. 



VII. AN ILLUSTRATIVE EXAMPLE 



In order to illustrate how the renormalization of the effective interaction works in a simple example, let us consider 
the dynamics of a point particle, diffusing in a rugged asymmetric harmonic oscillator: 

U(x) = h\x 2 + h%x + hs sm(wx), (70) 

with hi = 2, h-2 = 1, /13 = 1, tv = 4. The viscosity coefficient is set to 7 = 5 and inverse temperature to j3 = 5. Note 
that this potential has been chosen in such a way that the average value of x at thermal equilibrium is non- vanishing. 
The diffusion Monte Carlo algorithm used in our numerical simulations is presented in the appendix IS] The factor 



O, which appears in the L — 1 , L = 2 improvement terms was determined from the time interval At using Eq. ( 68 1 



The proportionality constant n in Eq. (69 1 was determined once and for all, by matching the result of (x(t)) of the 
unimproved (i.e. L = 0) simulations after 10 integration time steps with At = 0.01, with those of the RG-improved 
(i.e. L — 1, L — 2) MC simulations after a single elementary time step, with At' — 0.1. We found k = 0.35, with no 
appreciable difference between the L = 1 and L = 2 estimates. 

Let us now discuss the results of our simulations. We begin by analyzing the effects of accounting for the factor W(x) 
defined in Eq. (66 1, in numerical MC simulations. Fig. [5] shows the average position, once the system has attained 
thermal equilibrium, obtained by diffusion MC simulations with and without branching the walkers according to 
W(x). We recall that neglecting such a term is equivalent to simulating the dynamics in the Ito calculus, while 
the branching is expected to improve the time discretization to order At 2 . Indeed, our results show that, when one 
chooses small discretization steps, the two approaches are consistent with each other and yield the exact equilibrium 
average — which was computed directly from the Boltzmann distribution — . On the other hand, at large discretization 
steps, accounting for the factor W significantly improves the result. The same discussion can be trivially repeated 
in simulations in which the factor W(x) is interpreted as a re-weighting term, while the number of walkers is held 
constant. 

We now discuss the use of our effective theory to simulate the stochastic dynamics, using large time steps. Fig. [6] 
shows the time evolution of the average particle position at time t, computed using a small discretization time step 
— At = 0.01 — and a large discretization step — At' = 0.1 — . The two curves obtained in the original — i.e. tree- 
level — theory are compared with the results of the effective theory at order L = 1 and L = 2, which were obtained 
using an integration time step which was one order of magnitude larger, A t' = 0.1. 

The time evolution of the observable (x(t)), obtained in the tree-level theory using large integration time steps 
(squares) is inconsistent with the same quantity obtained using small time steps At = 0.01 (circles). This is expected, 
because for At = 0.1 the numerical simulations of the tree- level theory start to be affected by significant discretization 
errors — see Fig. [5]—. 
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Renormalization Group Improved MC 

Time evolution of the average position 




FIG. 6: The average position of the particle at time t, computed in the tree-level theory (At = 0.01 for L = 0), and in the 
effective theory (At = 0.1 for L = 1 and L — 2). The insert displays a part of the same curve, on a larger scale. Statistical 
errors are smaller than the symbols. The At — 0.1, L = 2 cannot be distinguished from the At — 0.01, L = curve. 

The results of simulations with large discretization time steps are significantly improved if one uses the effective 
theory, already at order L = 1 (diamonds). At order L — 2 the dynamics of the tree-level theory simulated at 
At = 0.01 is indistinguishable from the dynamics of the effective theory simulated with At' = 0.1 (triangles). These 
results show that the hard-mode dynamics in the short time range from 0.01 to 0.1 has been correctly taken into 
account by means of the renormalized effective interaction. As a consequence, the use of the effective theory allows 
to obtain very accurate predictions, using larger time steps. 

VIII. LONG TIME DYNAMICS OF MOLECULAR SYSTEMS 

The improvement of the MC algorithm based on our effective theory is expected to be most efficient when the gap 
between the slow and the fast modes is very large. In fact, in this regime, the slow-mode perturbation theory remains 
reliable even when one integrates out a large frequency shell, i.e. when &< 1. Hence, in this case, by RG-improvcmcnt 
it is possible to simulate the time evolution using elementary time steps At' which are significantly larger than the 
original elementary time step At, which would be used in the usual (unimproved) MC algorithm. 

A natural application of the RG-improved MC is the investigation of the long-time dynamics of macromolecules, 
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for which standard MD or MC algorithms can be extremely computationally expensive. Hence, it is interesting to 
address the question of what is the typical range of reliability of the slow-mode perturbation theory for a typical 
molecular interaction, at room temperature. To this end, let us consider the over-damped diffusion at temperature 
300 K of two molecules of mass to = 30 amu, interacting through a Van-Der-Waals potential: 



U(r) = 4e 



(?)"-(? 



(71) 



where e = 4 KJ/mol and a = 0.3 nm. A typical value for the viscosity coefficient for a molecule in its solvent (e.g. 
an amino acid in water) is 7 ~ 2 x 10 3 amu ps _1 .The typical time-steps used in the numerical integration of the 
Langevin Eq. (|2j) are of the order At ~ 10~ 3 — 10~ 2 ps. 



The tree-level effective interaction associated to the potential (71 1 is: 



Veffir) = 



4", 



24e 



1-2- 



8e a" 

~/3 r8 



156- 



42 



(72) 



This function and the corresponding L = 1 and L = 2 renormalized effective interactions are plotted in Fig. [7] for 
k = 1. This plot shows that, for a realistic set of parameters, the perturbative expansion remains reliable even when 
one integrates out a very large shell of modes, with b ~ 10 -2 . This fact suggests that the ultra-violet dynamics is 
essentially free brownian motion, while the long time dynamics is dominated by very low-frequency modes, and is 
driven by the force field. This fact has remarkable consequences on practical numerical simulations. It implies that 
by using the renormalized effective potential, it should be possible to adopt integration time steps which are about 
10 2 time larger than those required to simulate the dynamics in the original tree-level theory. 



IX. CONCLUSIONS 



In this work, we have presented a new approach to the problem of investigating the long-time out-of-equilibrium 
dynamics of multi-dimensional systems obeying Langevin dynamics. In the presence of decoupling of time scales, the 
methods based on the direct integration of the Langevin Eq. (MD) or on the time propagation of the Fokker-Planck 
probability density (MC) are usually inefficient, because a significant amount of computational time in invested to 
simulate uninteresting fast stochastic fluctuations. 

We have shown that the decoupling of time scales which limits MD and MC approach can in fact be exploited 
to perform analytically the average over the short-time stochastic fluctuations. After the integration over the fast 
modes has been performed, one obtains an effective theory which describes directly the relevant dynamics, with a 
lower time resolution. In such an effective theory, the effective action in the path integral receives corrections, which 
account for the ultra-violet physics which is cut-off. We have developed a rigorous scheme which allows to organize 
such corrections in term of a perturbative series in which the expansion parameters are the ratio between the soft 
frequency scales and the hard frequency scale b£l. Hence, sub- leading terms in the perturbative expansion come with 
higher inverse powers of the hard scales 6f2 and become irrelevant in the limit in which the decoupling of fast and 
slow modes is very large. 

The Feynmann diagrams which have to be calculated to obtain the corrections to any given order in this perturbation 
theory can be identified from their degree of slowness 

L(Feynman diagram) = N v — 1 + N T H — ^ (73) 

Diagrams with degree of slowness L generate corrections proportional to 1/(M7) L . In particular, we have found that 
the leading-order correction (i.e. L = 1) is proportional to the Laplacian of the effective potential V e ff. 

S>[X<] ~ A ^//(Mt)). (74) 

At the next-to-leading order, a term containing fourth-order derivatives appears: 

On the other hand, a space-dependent, tensor correction to the diffusion coefficient appears only as a higher-order 
effect (L = 5). It is important to stress the fact that, in the present approach, the ultraviolet cut-off f2 (or, equivalently, 
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FIG. 7: The tree-level, L= 1 and L = 1 renormalized effective potential 14// (r) for the Van-der-Waals interaction Eq. (71 1 
obtained integrating out the modes in the shell Sb = (bQ, fi) with SI = 27r/0.01ps and b — 0.01. 



the short time scale At) is kept finite at all stages. Upon taking the continuum limit At — » 0, all the correction terms 
in the effective theory vanish and one recovers the original theory, defined by the effective Schrodinger Eq. ([6| . 

The main usefulness of such an effective theory resides in the fact that it can be used to develop an improved MC 
approach, to compute the long-time evolution of the Fokker-Planck probability. The elementary time steps used in 
the RG improved MC algorithm are a factor 1 jb larger those of the MC algorithm for the underlying tree-level theory. 
Since the dynamics in the time range (At, 1 /b At) is averaged analytically, the RG improved MC algorithm avoids 
investing computational time in simulating the fast-mode dynamics associated to local Brownian motion. 

In the specific case of molecular interactions at room temperature, we have shown that the perturbative approach 
remains reliable even when one integrates large frequency shells, with b ~ 0.01. This feature suggests that, by using 
the effective theory, it is possible to simulate time intervals which can be up to a factor ~ 100 longer than in the usual 
MC approach. 

APPENDIX A: PROPAGATOR OF THE FAST MODES 

Here we derive free fast mode propagator G>(cj n , u> m ) appearing in the diagrams, using the standard source tech- 
nique. We first add a source term to Z> : 



z > -» Z°[v(" n )] 



Pis e 



-0 t Ei, 



x(tu n ) x(~-u) n ) -j- x (uj n )r}(— LJ n ) 



(Al) 



t 

2 Z^i^ies,, 



Then, we functionally differentiate twice with respect to the source: 



G°(w„,u; m ) = lim 



1 



lK)l(-"n) 



v^o (p t) 2 Sr)(-u n ) 8r](-w m ) 



(3 7 t ul 



(A2) 



(A3) 



Note that since the zero mode belongs to the slow modes part of the kinetic action, the kinetic operator for the fast 
modes is never singular and can be inverted without troubles. 



APPENDIX B: DIFFUSION MONTE CARLO ALGORITHM 



Our numerical study were performed using the following diffusion Monte Carlo algorithm: 
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1. A ensemble of N w = 18000 initial configurations {x\(t — 0), . . . , x^ w (t = 0)} was generated by sampling from a 
narrow Gaussian distributions of width a = 0.01, centered at the origin x — 0. Each of such positions represents 
the starting point of a walker. 

2. A new set of N w configurations was obtained by evolving the initial points for an elementary interval A t , 
according to the Langevin dynamics in the Ito calculus: 

xi(t + M) = xi(t)-— 4- U(xi(t)) + Atr], l = l,...,N w . (Bl) 

7 ax 

Atr/ represents the usual Brownian diffusion term, which was performed by sampling from a Gaussian of width 
a 2 = j- At, centered at the origin. 

3. For each walker, we generated a random number £ £ [—0.5, 0.5] and we made N c copies of the walker, where N c 
is the integer part of W(x(t + At)) + £. Hence, for N c = the walker was aborted, for N c = 1 the walker was 
left unchanged, while for iV c > 1 the walker gave raise to descendents, which then propagated independently 
from the progenitor. The integration time step At was chosen in such a way that the relative fluctuations in 
the population of walkers was only occasionally exceeding 10%. 

4. The steps 2-3 were iterated for many integration time steps. 

5. The quantity (x(t)) was obtained from the mean over the configurations of the walkers. The statistical error 
was estimated from the variance. 
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